宏基因组实战5. sourmash基于Kmer比较数据集
前情提要
如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章
测试数据
刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。
最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。
sourmash比较数据集
https://2017-cicese-metagenomics.readthedocs.io/en/latest/sourmash_compare.html
sourmash教程
主页：https://sourmash.readthedocs.io/en/latest/
Mash: fast genome and metagenome distance estimation using MinHash. Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, Phillippy AM. Genome Biol. 2016 Jun 20;17(1):132. doi: 10.1186/s13059-016-0997-x.
sourmash是一款超快、轻量级核酸搜索和比对软件，方法叫MinHash。这这一个命令行使用的python包，可以从DNA序列中快速分析kmer，并进行样品比较和绘图。目的是快速且准确的比较大数据集。在2016年6月20日发表，目前已经被引用61次。
目的
比对reads至拼装结果
比较数据集并构建聚类图
采用k-mer Jaccard distance进行样品比较
什么是k-mer
k-mer就是长度为k的DNA序列
ATTG - a 4-mer
ATGGAC - a 6-mer比如，一个长度为16的序列可以提取11个长度为6的 k-mers
AGGATGAGACAGATAG
AGGATG
GGATGA
GATGAG
ATGAGA
TGAGAC
GAGACA
AGACAG
GACAGA
ACAGAT
CAGATA
AGATAG需要记住的概念：
不同物种的k-mer是很不同的
长k-mer具有很强的物种特异性
K-mers与组装图
三大基因组拼接方法之一，就是将基因组打断成k-mer再拼接，通过k-mer步移实现拼接
为什么采用k-mers而不是全长序列拼接
计算机喜欢k-mer，因为匹配准确快速。
长k-mers存在物种特异性
图1. 以k=40为例，kmer很容易按属水平分开细菌
采用k-mers比较样品
安装sourmash
如果你进行过之前几节课的分析，只需执行下面代码的最后一条即可
# 设置工作目录，这是我的目录，学习时请修改为你工作的目录
pwd=~/test/metagenome17
cd ${pwd}
# 依赖关系，之前安装成功可跳过
sudo apt-get -y update && \
sudo apt-get install -y python3.5-dev python3.5-venv make \
libc6-dev g++ zlib1g-dev
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer
# 安装程序包
pip install -U https://github.com/dib-lab/sourmash/archive/master.zip产生Illumina reads的signature
mkdir work
cd work
curl -L -o SRR1976948.abundtrim.subset.pe.fq.gz https://osf.io/k3sq7/download
# 此文件很大，如果继续以之前的分析，可以执行下行链接至此目录
# ln ../data/SRR1976948.abundtrim.subset.pe.fq.gz ./
curl -L -o SRR1976948.megahit.abundtrim.subset.pe.assembly.fa https://osf.io/dxme4/download
curl -L -o SRR1976948.spades.abundtrim.subset.pe.assembly.fa https://osf.io/zmekx/download图2. sourmash工作流程
mkdir ../sourmash
cd ../sourmash
# 计算过滤序列的k51特征
sourmash compute -k51 --scaled 10000 ../work/SRR1976948.abundtrim.subset.pe.fq.gz -o SRR1976948.reads.scaled10k.k51.sig比较序列与组装结果
sourmash compute -k51 --scaled 10000 ../work/SRR1976948.spades.abundtrim.subset.pe.assembly.fa -o SRR1976948.spades.scaled10k.k51.sig
sourmash compute -k51 --scaled 10000 ../work/SRR1976948.megahit.abundtrim.subset.pe.assembly.fa -o SRR1976948.megahit.scaled10k.k51.sig图3. 评估污染比例
sourmash search SRR1976948.reads.scaled10k.k51.sig SRR1976948.megahit.scaled10k.k51.sig --containment
sourmash search SRR1976948.reads.scaled10k.k51.sig SRR1976948.spades.scaled10k.k51.sig --containment结果如下：
loaded query: SRR1976948.abundtrim.subset.pe... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
48.7% SRR1976948.megahit.abundtrim.subset.pe.assembly.fa
loaded query: SRR1976948.abundtrim.subset.pe... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
47.5% SRR1976948.spades.abundtrim.subset.pe.assembly.fa为什么只有不到一半的Kmer在拼接结果中？
我们看看拼接结果中有多少kmer在原始序列中
sourmash search SRR1976948.megahit.scaled10k.k51.sig SRR1976948.reads.scaled10k.k51.sig --containment
sourmash search SRR1976948.spades.scaled10k.k51.sig SRR1976948.reads.scaled10k.k51.sig --containment结果如下：
loaded query: SRR1976948.megahit.abundtrim.s... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
99.8% SRR1976948.abundtrim.subset.pe.fq.gz
loaded query: SRR1976948.spades.abundtrim.su... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
99.9% SRR1976948.abundtrim.subset.pe.fq.gz基本都全部可以找到。这是因为原始序列中包含大量illumina测序的错误，而在拼接中被丢弃或校正。
特征比较
为了更深刻理解，我们采用osf下载更多数据集比较
pip install osfclient
osf -p ay94c fetch osfstorage/signatures/SRR1977249.megahit.scaled10k.k51.sig SRR1977249.megahit.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977249.reads.scaled10k.k51.sig SRR1977249.reads.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977249.spades.scaled10k.k51.sig SRR1977249.spades.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.megahit.scaled10k.k51.sig SRR1977296.megahit.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.reads.scaled10k.k51.sig SRR1977296.reads.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.spades.scaled10k.k51.sig SRR1977296.spades.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/subset_assembly.megahit.scaled10k.k51.sig subset_assembly.megahit.scaled10k.k51.sig图4. 样品间比较原理
# 比较所有signature文件
sourmash compare *sig -o Hu_metagenomes
# 比较结果绘图
sourmash plot --labels Hu_metagenomes生成Hu_metagenomes.dendro.png和Hu_metagenomes.matrix.png两个文件，分别是树形图，还有热图+树形图，可以用Filezilla下载查看。
图5. 样品间kmer聚类结果
此软件还可以做什么？
更多的kmer研究
物种分类
功能分析
想学习看该软件的官方帮助文档吧。https://sourmash.readthedocs.io/en/latest/
猜你喜欢
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外七十多位PI，七百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读